Chapter 12 High-level functions

In most situations you will have rasters that are larger than what you need. You might be interested in the landcover of your study area but the data covers all of Europe. Conversely, you might also need to combine different rasters if the scale of your analysis is larger than the data you are using. Often this is the case when larger rasters are split in smaller areas to reduce their size.

In the first case you would use crop(). Provided with a raster and an extend, crop() will return only the part of the raster that lies within the extend. An extend usually is a named vector with four elements: xmin, xmax, ymin and ymax. It might for example be the bounding box you derived from the vector data you were analyzing in sf. You can create this raster yourself if you know, or can derive, the relevant numbers. An alternative is the drawExtent() function from the raster package. This lets you click on an already plotted raster an returns the extend of the thus created rectangle. Try the code below in your own R instance to see what I mean.

plot(dem1)
extent = raster::drawExtent()
dem1crop = terra::crop(x =dem1, y=extent)

Nice! Now have reduced the raster to the area we are actually interested in. We can compute the mean hight or the sum of all hights with global().

global(dem1crop, "sum")

##                    sum
## DGM25_2905530 100724.5

global(dem1crop, "mean")

##                   mean
## DGM25_2905530 331.3307

We can also expand a raster with the extend() function. That means we add rows and/ or columns with NA values. The other way around we can also trim them, removing all cells with NAs in the margins.

dem1.exp = extend(dem1, y = c(100,100))
plot(dem1.exp)

Since all the NAs are just white, it’s hard to see what we actually did here. Let’s give them some random value so we can see the cells we added.

na_cells = which(is.na(values(dem1.exp)))
values(dem1.exp)[na_cells] = 300
plot(dem1.exp)

We can now clearly see the added cells in brown. With trim() we can remove them again.

na_cells = which(values(dem1.exp) == 300)
values(dem1.exp)[na_cells] = NA
dem1.trm = trim(dem1.exp)
plot(dem1.trm)

As mentioned before, sometimes you want to combine rasters. For example our two DEMs are both part of a larger DEM that covers all of the German federal state Rhineland-Palatinate. We can combine the two rasters with merge()

dem_merge = merge(dem1, dem2)
plot(dem_merge)

Whenever you want to use multiple raster in an operation, it is important that they have the same CRS and cell sizes. You can increase cell sizes with aggregate() and reduce cell sizes with disaggregate(). In the first case, several cells are combined into one. The number of cells is set with the fact argument. When you supply one number it is used for both directions. So if you call aggregate(x, fact = 2) and raster has the cell size 10m by 10m, the new raster will be 20m by 20m. You can also provide fact with a vector that specifies the aggregation factor in each direction. The third argument is fun which specifies how to calculate the new cell value. If our new cell consists of four old cells (i.e., we used an aggregation factor of 2) we somehow need to combine those four numbers into one. Common options here are the mean, the median, the modal, the maximum, the minimum, or the sum.

dem1_agg_2_mean = aggregate(dem1, fact=2, fun = "mean")
dem1_agg_3_mean = aggregate(dem1, fact=4, fun = "mean")
dem1_agg_4_mean = aggregate(dem1, fact=6, fun = "mean")
dem1_agg_5_mean = aggregate(dem1, fact=8, fun = "mean")
par(mfrow = c(2,2))
plot(dem1_agg_2_mean)
plot(dem1_agg_3_mean)
plot(dem1_agg_4_mean)
plot(dem1_agg_5_mean)

par(mfrow = c(2,2))
dem1_agg_mean = aggregate(dem1, fact = c(1,6), fun = "mean")
dem1_agg_min = aggregate(dem1, fact = 4, fun = "min")
dem1_agg_max = aggregate(dem1, fact = 4, fun = "max")
dem1_agg_sum = aggregate(dem1, fact = 4, fun = "sum")
plot(dem1_agg_mean,  main = "directed aggregation")
plot(dem1_agg_min,   main = "minimum")
plot(dem1_agg_max,   main = "maximum")
plot(dem1_agg_sum , main = "sum")

When we want to go in the other direction, we disaggregate raster cells with the disagg() function. Sadly, we cannot magically determine what values are the field truth for the new, smaller cells. Instead each new cell gets the same value as it’s parent.

dem1_disagg = disagg(dem1, fact = 2)
plot(dem1_disagg)

12.1 Extracting raster values with vectors

We often encounter situations where we want to extract the values of a raster that coincide with vector data we have. For example, we might have taken soil samples along a mountain and now we want to now the altitude of the samples or we are interested in the land cover of a river catchment.

To demonstrate this, we will need some vector data sets. For the point data, we can simulate points within the bounding box of the DEMs we used before.

# - extract extend of DEM
dem.ext <- ext(dem1)
# - convert extend to bounding box usable by sf
dem.bb <- st_bbox(dem.ext)
# - create 10 random points within bounding box
random_points <- st_sample(x = st_as_sfc(dem.bb), 10)
random_points <- st_as_sf(random_points, crs = 25832)
#- interactive map
tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(dem1) + tm_raster() + tm_shape(random_points) + tm_dots(size = 1)

Now we have ten points that are randomly distributed across the map. If your repeat this with you own computer your ten points will be at other positions. Now that we have, we can extract the raster values, i.e., the altitude, at the locations of our points. This extraction is done with extract(). The first argument to extract() is the raster raster from which we want to extract values. The second argument is the location of points, at which we want to extract the values.

ex.points <- extract(dem1, random_points)

We get a non-spatial data.frame back that holds the row number of points and the elevation at these points. We can add the altitude values back to the points like this:

random_points$altitude <- ex.points$DGM25_2905530
tm_shape(random_points) + tm_dots(col = "altitude")

Now we want to look at polygons instead of points. First, we need some polygons. Here, we use river catchments from Rhineland-Palatinate.

catchments <- st_read("data/catchments.gpkg")

## Reading layer `catchments' from data source `C:\Users\jonat\Documents\001_Uni\002_teaching\001_GIS\gisbook2\data\catchments.gpkg' using driver `GPKG'
## Simple feature collection with 2754 features and 55 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4042743 ymin: 2874149 xmax: 4212789 ymax: 3094667
## Projected CRS: ETRS89-extended / LAEA Europe

Let’s have a look at the data.

tm_shape(catchments) + 
        tm_polygons(alpha = .5) + 
        tm_shape(dem_merge) + 
        tm_raster()

As we can see, we are missing elevation data for most of the catchments. Thus, we will first remove all catchments for which we are lacking data.

dem.ext <- ext(dem_merge)
# - convert extend to bounding box usable by sf
dem.bb <- st_bbox(dem.ext) |> st_as_sfc() |> st_as_sf(crs = 25832)
# - transform catchments to same bb as the DEMs
catchments <- st_transform(catchments, crs = 25832)
# - subset catchments to only those that are completely within the DEM 
catchments2  <- catchments[dem.bb, , op = st_within]

Let’s have a look at the result.

tm_shape(catchments2) + 
        tm_polygons(alpha = .5) + 
        tm_shape(dem_merge) + 
        tm_raster()

Next, we extract the raster values from the cells within the different catchment polygons with the extract() function. There is a function with the same name in the raster package. To make sure that we are using the function from terra, we should therefore call the package explicitly by using terra:: in front of the extract() function.

dem_ext <- terra::extract(dem_merge, catchments2)
head(dem_ext)

##   ID DGM25_2905530
## 1  1       458.403
## 2  1       458.774
## 3  1       459.436
## 4  1       460.564
## 5  1       461.968
## 6  1       463.475

The result is a data.frame with two columns: ID and DGM25_2905530. The second cell is named after the raster from which we extracted the values. The ID column assigns the raster values to the polygons we used for the extraction.

unique(dem_ext$ID)

## [1] 1 2 3 4 5 6 7

As we can see, the ID values go from one to seven and we have

nrow(catchments2)

## [1] 7

seven catchment polygons. So all rows in dem_ext in which the ID variable is 1 contain raster cell values from cells that are within the polygon in the first row on catchments2.

We could now compute the mean elevation within each catchment and assign those values as a variable to catchments2. Here, I will show two different ways to summarize one variable by another. In our case, the elevation (DGM25_2905530) by the catchment (ID). First, we can use group_by() and summarize() from the dplyr package. The first function defines the variable that groups the rows and the second computes some function of a variable.

dem_ext |> 
        group_by(ID) |> 
        summarize(mean_elev = mean(DGM25_2905530))

## # A tibble: 7 × 2
##      ID mean_elev
##   <dbl>     <dbl>
## 1     1      392.
## 2     2      283.
## 3     3      377.
## 4     4      327.
## 5     5      278.
## 6     6      385.
## 7     7      364.

Alternatively, we can use the data.table package. Here, we need to turn dem_ext into a data.table object with setDT(). setDT() is a mutating function. That means we do not need to assign the result to an object but the object within the function is changed (mutated) directly.
Data tables can be subset or modified following the basic scheme dt[i,j,by], we i is a logical expression subsetting the rows, j is an expression subsetting or modifying columns and by provides a grouping structure for j.

library(data.table)
setDT(dem_ext)
(dem_ext2 <- dem_ext[, mean(DGM25_2905530), by = "ID"])

##    ID       V1
## 1:  1 391.8343
## 2:  2 282.9041
## 3:  3 377.3319
## 4:  4 326.7337
## 5:  5 277.8304
## 6:  6 385.4874
## 7:  7 363.8802

Either way we get the same results, which we can now add back to the catchments.

catchments2$elev_new <- dem_ext2$V1
tm_shape(catchments2) + tm_polygons(col = "elev_new")

12.2 Masking

Sometimes we might want to fill holes in rasters, i.e., areas that only contain NAs or some other non-informative value, with the values from a second raster. This operation is called masking. Here we will work through a short example, where there are randomly scatted NAs throughout the raster. This might happen if small but thick clouds prevented our satellite to accurately measure conditions on the ground.

# Example1 
dem1_na = dem1 
values(dem1_na)[sample(1:ncell(dem1), size = 100)] = NA
#to increase the visibility of the masked cells we will create a separate mask layer with extreme values. 
dem1_mask = dem1
values(dem1_mask) = 1000
dem1_nax = cover(x = dem1_na, y = dem1_mask)
plot(dem1_nax)

The lapp() function is used for functional programming. It is similar to apply()-family functions or map-family functions from the purrr package. It applies a function to each layer of a raster Here for example we increase the height of a third of the cells by 10m.

dem1_2 = dem1_3 = dem1
values(dem1_2) = rbinom(ncell(dem1),1,0.33)
values(dem1_3) = 10
dem1_all = c(dem1, dem1_2, dem1_3)
dem_lapp = lapp(dem1_all, function(x,y,z) x + y * z)
plot(dem_lapp)

We can classify the cells with the classify() function.

m <- c(0, 100, 1,
       100, 250, 2,
       250, 400, 3,
       400, 600, 4)
dem1_cut = classify(dem1, rcl = m)
plot(dem1_cut)

focal() can be used to replace the value of a focal cell by some function of its neighbours. Which neighbors and what function to use can be chosen as arguments.

mw = focal(dem1, w=3)
plot(mw)